####========================== Libraries ==============================####

for (library_name in c("ggthemes", "RColorBrewer", "gtable", "stargazer",
                       "usmap", "ggplot2", "gridExtra", "maps", "lubridate", "stringi", "stringr", "RCurl",
                       "sandwich", "lmtest", "tidyverse", "readxl", "xtable", "gender", "wru",
                       "reshape2", "reldist", "shades", "colorspace", "ggstance", "cowplot", "haven",
                       "rmarkdown", "usmap", "usmapdata", "ggpattern", "xlsx")) {
  library_name
  if (!require(library_name,character.only=TRUE)) {
    install.packages(library_name, character.only=TRUE)
  }
  library(library_name, character.only=TRUE)
}


####============================= Setup ================================####
# these need to be changed to your own directory layout
tables_path <- "/data"    
output_path <- "/output"   

####============================= Figure 4 Maps ================================####
data <- read_dta(file.path(tables_path, "state_year_long_avgs.dta"))

all_states <- usmap::statepop %>%
  rename(state = abbr) %>%
  subset(select = state)

data <- data %>%
  filter(state != "USA" & state != "USA_balanced") %>%
  merge(all_states, by = "state", all.x = TRUE, all.y = TRUE) %>%
  mutate(y0_cuts = factor(case_when(re_inv_y0_uw_mean < -.01 ~ "< -.01",
                                    between(re_inv_y0_uw_mean,-.01,0) ~ "-.01 - 0",
                                    between(re_inv_y0_uw_mean,0,.01) ~ "0 - .01",
                                    between(re_inv_y0_uw_mean,.01,.05) ~ ".01 - .05",
                                    re_inv_y0_uw_mean > .05 ~ "> .05"),
                          levels=c("> .05",".01 - .05","0 - .01",
                                   "-.01 - 0","< -.01")),
         y1_cuts = factor(case_when(re_inv_y1_uw_mean < -.01 ~ "< -.01",
                                    between(re_inv_y1_uw_mean,-.01,0) ~ "-.01 - 0",
                                    between(re_inv_y1_uw_mean,0,.01) ~ "0 - .01",
                                    between(re_inv_y1_uw_mean,.01,.05) ~ ".01 - .05",
                                    re_inv_y1_uw_mean > .05 ~ "> .05"),
                          levels=c("> .05",".01 - .05","0 - .01",
                                   "-.01 - 0","< -.01")))
         
         plot_usmap(data = data, regions = "states", value = "y0_cuts", linewidth = 0.2) +
           scale_fill_manual(values = c("< -.01" = "#2a9cf4",
                                        "-.01 - 0" = "#ccefff",
                                        "0 - .01" = "#faf591",
                                        ".01 - .05" = "#ffc334",
                                        "> .05" = "#d15000"),
                             name = "UD (Y*=0) Mean") +
           theme(legend.position = "right") +
           geom_polygon_pattern(
             data = us_map(regions="states") %>%
               rename(state = abbr) %>%
               right_join(data, by = "state") %>%
               filter(nodata == 1),
             aes(x,y,group=group),
             color = NA,
             linewidth = 0.2,
             fill = 'white', 
             pattern_spacing = 0.05, 
             pattern_density = 0.3, 
             pattern_fill = "black")
         ggsave(file.path(output_path, "y0_mean_r.pdf"), device = "pdf", width = 8, height = 4.5, dpi = 300)
         
         plot_usmap(data = data, regions = "states", value = "y1_cuts", linewidth = 0.2) +
           scale_fill_manual(values = c("< -.01" = "#2a9cf4",
                                        "-.01 - 0" = "#ccefff",
                                        "0 - .01" = "#faf591",
                                        ".01 - .05" = "#ffc334",
                                        "> .05" = "#d15000"),
                             name = "UD (Y*=1) Mean") +
           theme(legend.position = "right") +
           geom_polygon_pattern(
             data = us_map(regions="states") %>%
               rename(state = abbr) %>%
               right_join(data, by = "state") %>%
               filter(nodata == 1),
             aes(x,y,group=group),
             color = NA,
             linewidth = 0.2,
             fill = 'white', 
             pattern_spacing = 0.05, 
             pattern_density = 0.3, 
             pattern_fill = "black")
         ggsave(file.path(output_path, "y1_mean_r.pdf"), device = "pdf", width = 8, height = 4.5, dpi = 300)